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FFT-R - A FAST FOURIER TRANSFORM SUBROUTINE 
FOR REAL VALUED FUNCTIONS 


* 


DECUS Program Library Write-up 


DECUS No. 8-143 


1. ABSTRACT 

The Fast Fourier Transformation enables computation 
of the power spectrum of a time series in a minimum of 
time. Specifically, it reduces the number of computations 
required to calculate the Discrete Fourier Transformation 

5 i = ^r2L^k wJk Cw--c 2T ' i/ 'V=V=T) 

J N iTo k 

of a series of N equally time spaced samples Xo, Xj> ,. 

Xn- 1 where N is a power of 2(N=2 n ). in fact, for 1024 
time samples, computation time is reduced by over 99%. 

FFTS-R (for Fast Fourier Transformation Subroutine) 
will transform up to 2048 real points. It is written as 
a subroutine, and is I/O independent. The user must tailor 
his own input-output procedure to his particular environment. 

2. REQUIREMENTS 

2.1 Hardware 

A 4K PDF-8 with Extended Arithmetic Element Type 182 
or a PDP-8/I with EAE Type KE-8/l option is the 
minimum necessary hardware. 

2.2 Storage 

FFTS requires locations 3 to 7, 20 to 107, and 400 to 
2401+N, where N is the (octal) number of points being 
transformed. 

3. LOADING PROCEDURE 

Make sure the BIN Loader is in core. If not, load it. 
Put 7777 in the SR. Press Load Address. Place FFTS on the 
reader and turn the reader on. Press start, and FFTS will 
load. Then load the user's program the same way as above 
and start it. 
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4. USAGE 


4.1 Calling Sequences 


FFTS enables the user to take either the Fast Fourier 
Transform, (FFT) or its inverse (IFFT) of a real valued 
time series. The subroutine FFT, which begins at 0400, 
calculates the FFT. Register DOFFT (location 0060) points 
to FFT, so a JMS I DOFFT (=4460) will call FFT. The sub¬ 
routine IFFT beginning at 0076 takes the inverse FFT. 

Since location DOIFFT (normally 0061) points to IFFT, IFFT 
can be executed simply by writing JMS I DOIFFT (=4461). 
Both FFT and IFFT assume that the real data to be handled 


has already been stored in memory (see section 5). 

After the operation is complete, the results will be 
stored in memory in bit inverted order (see section 5.1). 
For FFT, the results are the complex co-efficients S . 
(with the appropriate scale factors, as described in J 
section 5.2) given by the equation in section 1 (j=0, 
1,.#., N-l). For IFFT the results consist of a time 
sequence Xj (j=0, 1,....,N-1). 


NOTE s THE REMARKS IN THE FOI ,LOWING SECTIONS APPLY TO 
IFFT AS WELL AS FFT. 


4.2 Execution Times 


The following is a table of execution times for the 
subroutine. 


Number of points transformed Time (Seconds) 


2046 

1024 

512 

256 

128 

64 


4.95 

2.20 

.963 

.417 

.177 

.074 


5. DETAILS OF STORAGE 
5.1 Data Storage 

A JMS I DOFFT causes a real time series to be Fourier 
transformed. That series is stored in memory. More 
explicitly, the data is stored sequentially after 
location XRTAB (=2400). For example, the storage scheme 
for a N=8 point transformation would be as follows 
(Xi is the i^h time sample)s 
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XRTAB, 


*2400 

X 0 

*1 

X 2 

X 3 

x 4 

X 5 

x* 


On exit the results of the transformation will be in 
core, only half of them, however, will actually be 
present. This is because the program makes use of a 
Hermitian symmetry in the frequency domain to save both 
time and storage space. The symmetry is as follows: 

If the time sequence x 0 , X 1 ,...X N _ 1 is real valued, 
then the pair of Fourier co-efficients s .• and S M * 

obey a complex conjugate symmetry. That J is s i = J S N . 
where '*' denotes taking the complex conjugate. N ” J 

So due to this symmetry, only one half of the co-cffi- 
cients need actually be computed, since either half is 
derivable from the other. Hence FFTS computes only Sq 
through S N / 2 / introducing a time and space saving factor 
of 2 , yet sacrificing no information. 

After execution the set of complex numbers Tsq,.. •/S N /oV 
are to be found in memory. Unlike the original data^' 
set they are not in sequential order, but 

rather in something called bit inverted order. Bit in¬ 
version means simply the process of re-ordering the bits in 
a binary number. For instance, the binary number 001 
bit inverted is just 100 (=4). Thus to locate the 
Fourier co-efficients s i (j<N/2), write j as a binary 
number of n-l(=log 2 (N) J bits, bit invert j, and look 

2 

in that position. For example, to locate S 2 in memory 
for an 8 point transformation (N=8, n=3, n-I=2 ) write 
2 as a binary number of n-l=2 bits, 2 10 =10 2 . Then 
reverse the order of these bits giving 01=1 10 . This 
means that S 2 is stored in position 1. Physically, 
then, S 2 is to be found in location 2400 + 2 (1). The 
reason that bit inverted j (= 1 ) is multiplied by 2 is that 
each Sj is complex, so two locations are required to 
store it - one for the real part, the other for the 
imaginary one. FFTS adopts the following format: the 
imaginary part of a number is stored in the register 
immediately following the real part. As a specific 


MEANS REAL PART 
MEANS IMAGINARY PART 


point transform is ■ 

written out 


*2400 


XRTAB, 

RE(S 0 ) 

/ RE 


• IM (Sq) 

re(s 2 ) 
im(s 2 ) 
RE(si) 
IMCSjl) 

re(s 3 ) 

IM(S 3 ) 

/ IM 


RE(S 4 ) 

im(s 4 ) 
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s N/2 (here S 4 ) always is placed last. If J denotes bit 
inverted j, then a general formula for locating the real and 
imaginary parts of Sj is (LOC ( ) denotes location of): 

LOC (RE (S j) ) =2400+2 J* 

LOC (IM(S j) )s2401 +*2j = LOC (RE(Sj)) + l 

where j is written out as a binary number of n -1 
bits. A subroutine INVRT (location 1036) has been 
provided to do bit inversion. It can be called by 
a JMS I INVERT (INVERT=55). See section 8 . 


In addition^ a subroutine SORTX has been included 
which sorts the co-efficients and leaves them in 
sequential order. It can be called by a JMS I SORT 
(SORT=54). If SORTX were called after an 8 point 
transform had been completed, the data buffer would 
look like this: 

*2400 

XRTAB, RE (So) 

IM(S 0 ) 

RE (Si) 

IM(S2) 

RE(S2) 

im(s 2 ) 

re(s 3 ) 

IK (S 3 ) 

RE (S 4 ) 

IM(S 4 ) 

The reason that the co-efficients are not automatically 
sorted is that time can be saved by outputing from bit 
inverted order, and this possibility should be allowed 
for. 


5.2 Data Scaling 


All calculations in FFTS are done with single precision 
fixed point signed binary fractions. The binary point 
is located between bit 0 and bit 1 , leaving an 11 bit 
signed mantissa. Bit 0 is used as a sign bit. Negative 
numbers are formed by taking the two's complement of 
the positive binary fraction. So all inputs must be 
scaled in magnitude to less than one. The outputs are 
also formatted as above. There is also a more subtle 
scale factor involved, in order to utilize the 
maximum number of bits in the transformation it is some¬ 
times necessary to divide by 2 in a computation. As 
a result of this a pseudo floating point format has 
been adopted in which a variable scale factor (or 
exponent) is imposed on all the Fourier co-efficients. 
This scale factor or pseudo exponent is found in regis¬ 
ter SCALE (= 66 ) after each transform has been completed. 
The numbers stored in memory are the Fourier co¬ 
efficients multiplied by 2 raised to the contents of 
SCALE. So to retrieve the co-efficients themselves, 
merely shift each number C(SCALE) places right. If 
any further computations are to be done, better 
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accuracy will be obtained by retaining the pseudo 
exponent and leaving the co-efficients in "normalized 
form." in the case of the inverse transform, the 
desired results (here time samples) are the numbers 
stored in memory times 2l^i-c(SCALE)).* 


6 . RESTRICTIONS 

6.1 Program Initialization 

Because FFT is a subroutine certain registers must 
be primed before the first entry in order to insure 
proper operation. Specifically, register M (location 
0020 Jmust contain the number of points being trans¬ 
formed (in octal, of course) and register MU (location 
0021) must contain the power of two which M is, that 
^ is, C (M)=2f C (MU) .C (MU) must be at least 3 and no more 

than 13g, due to memory limitations. 

6.2 Input Restrictions 

So as to prevent overflow of the single precision 
storage, it is absolutely necessary that all data be 
less than 1 in magnitude , subject to the format 
described in section 5.2. (The binary point is to 
the right of bit 0). 


7. METHODS 


7.1 Algorithm 

FFTS uses the algorithm discovered by Cooley and Tukey 
for the rapid computation of a spectrum. This al¬ 
gorithm, called the Fast Fourier Transformation (or 
FFT), permits transformation of N (which must be an 
integer power of 2) equally spaced time samples in a 
time proportional to Nlog 2 N, whereas previous methods 
required times proportional to N . This gives a 
reduction of l-log 2 N/lJ. For Nr 1024, this is over 99%, 
essence, the algorithm makes use of the fact that 


yy k - yjCk"'** M) 


* The inverse is defined here to be 

J k- c 

(s j are real), without the 1/lsi scale factor. 


In 
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(where W=e ' ) to reduce the number of multiplica¬ 

tions necessary for a transformation. A complete 
description and proof of the algorithm used and its 
implementation can be found in an article by James Rothman 
which appears in DECUSCOFE, Volume 7, Number 3. 


DETAILS OF OPERATION 

The following is a list of useful subroutines and 
their operations; (values of the symbols may be found 
in the symbol table included in this document.) 


Name 

Call 


Functions 

FFT 

JMS I 

DOFFT 

Takes the Fourier Transforma¬ 
tion of the data buffer. 

Results in bit reversed order. 

IFFT 

JMS I 

DOIFFT 

Takes the Inverse Fourier 
Transformation of the data 
buffer. Results in bit 
reversed order. 

SORTX 

JMS I SORT 

c+JJsk int/ertT 

Sort the data buffer so that 
it is in normal sequence. 

TRIGET 

JMS I 

GETRIG 

Fetches sine and cosine values. 
Specifically, if the AC=K on 
entry, the values of sin 
(2irKA) and cos (2|TK/kl) are 
fetched from an internal trig 
table. K must be <or=N/2. A 
register COSINE contains the 
cosine value and the AC con¬ 
tains the sine value on exit. 

INVRT 

JMS I 

INVERT 

Number in AC is bit reversed 
and the result is in the AC 
on exit. 

MULTIP 

JMS I 

MULT 

Rounded single precision signed 
multiply. Uses EAE. AC= 
multiplier. C(Call address + 1) 
=address of multiplicand. 

Result in AC on exit. 


SYMBOL TABLE 
A symbol table follows 






symbol 

table 

SYMBOL 

table 

SYMBOL TABLE 

adder 

2053 

bult IP 

1020 

XRT AR 2420 

aqdr 

1134 

MUY 

7425 

XSGN 3074 

addwos 

1156 

N 

3023 

XTRACT 1242 

AUDXTR 

3760 

NMI 

7411 

XTRADO 3063 

ADDl 

1173 

NOROT 

3571 


Al)D2 

3041 

notnor 

1171 


4003 

3777 

NOVER4 

0022 


ADJSGM 

3575 

l\|0 4 M * * 

1132 


A 1 

3052 

NU 

3024 


AH 

3051 

P 

3032 


AHG2 

1017 

PI 

0030 


ASR 

7415 

PH 

0027 


B1 

3050 

Q 

3031 


HIGeMU 

0013 

Q1 

0026 


Bill 

3037 

OH 

3025 


bilr 

0036 

QUAD1 

1110 


binml I 

0035 

0UAD2 

1072 


binmlr 

0034 

RBUILD 

3 7 2« 


BH 

3047 

HtCHK 

3726 


BUILD 

0551 

RESETC 

0725 


C 

0040 

RLVERS 

0713 


CAM 

7621 

S 

0006 


CC I A 

3136 

SCA 

7441 


chkp t 

3524 

scale 

3066 


CiMOP 

3127 

SCL 

7423 


CNOTS 

3722 

SETC 

3547 


cos I ml 

30 4 4 

SGNADJ 

3075 


DATAHI 

6432 

SGNX 

1314 


DOFF T 

0060 

SHFCHK 

0070 


DOIFFT 

0061 

shflag 

3067 


L) V I 

7407 

shfti 

1077 


F 

0027 

SHFT2 

1114 


FFT 

0420 

SHF T 3 

1125 


flip 

1044 

SHIFCT 

3567 


FLIPCT 

1060 

SHIFT1 

0071 


GtTRIG 

3057 

SHIFT2 

0072 


61 

0046 

SHIFT3 

3073 


GR 

0045 

SHL 

7413 


IFFT 

3076 

SIGN 

1035 


INDEX 

1133 

SINE 

3043 


INVERT 

3055 

SlNLOC 

0064 


INVRT 

1036 

SINRET 

1122 


K 

3033 

sintab 

1375 


L 

30 35 

SORT 

0054 


LOOP1 

04 50 

SORT X 

3737 


LSR 

7417 

SWAPEO 

0751 


M 

0020 

SYNTH 

3062 


MAXNU 

0023 

s yn t ht 

1174 


MNOVR2 

3024 

tlmpr 

0042 


mqa 

7521 

TRIGET 

1061 


mol 

7421 

RORD 

1056 


MU 

0021 

HORDP 

1057 


mult 

0056 

XHLOC 

3065 
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ADDENDUM TO 8-143 and 8-144 


The program was structured so to make the change of eliminating 
the EAE requirement with a minimum of effort. 

All that need be done is replace each EAE instruction with a 
subroutine that performs the given operation using a pseudo 
multiplier-quotient. For this purpose the EAE simulator may be 
used. This does not allow certain microcodes, and where these 
occur in the FFT program, they can be separated into groups of 
EAE instructions, all of which together perform the designated 
function. 

For example CLA MQL MUY (microcoed) could become the three 
instructions: 

CLA 
MQL 
MQA. 
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CORRECTION TO DECUS NO. 8-143 AND 8-144 



ORIGINAL 


CORRECTED 

CHANGE 


*1000 


*1000 



MULT IP, 

0 

• 


MULTIP, 0 

• 




• 

• 

RAL 


• 

• 

io/-f RAR 

'Joto 

* 


DCA SIGN 


DCA SIGN 

f 



MUY 


MUY 



ARG2, 

HLT 

ARG2, 

HLT 




SHL 


SHL 




0 


0 




DCA ARG2 


DCA ARG2 




SHL 


Lot* TAD SIGN 


* 


0 


SHL 

^i3 

* 


MQL 


?r jZf 


* 


TAD SIGN 


ic TAD ARG2 

1 ?lj 

* 


CLL RAR 


V) SPA 

Qsio 

* 


TAD ARG2 


it CLA CLL CMA RAR 

yzso * 


MQA 


2 * NOP 


* 


SZL 


SZL 




CMA IAC 


CMA IAC 




JMP 1 MULTIP 


JMP 1 MULTIP 


SIGN, 

0 

SIGN, 

0 




The error was in the way in which rounding was accomplished. This fix was tested by performing a 
DOFFT, SORT, DOIFFT, SORT sequence on a 512 point real valued time series with 8-144 and then 
summing the absolute value of the imaginary residuals. The fix above reduced the sum by 40 percent. 





COBfiECTION TO DECUS NO. 8-143 


by 
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- 1 - 

The subroutine DECUS NO. 8-143 (may be DECUS 
NO, 8-144 also) does not consider an overflow which 
may occur during a complex multiplication: 

(a+ib).(cos O+i sin 0) * (a.cos 0-b sin 0)+i(a.sin©+b cosO) 

* c+id 

For 7F< ]a| + |b| < 2, an overflow is apt to occur 
for some values of 0. For example, if a = 0,9, b = 0,8, 

0 =Ti/4, then c ■ 0,07, d * 1.20, which causes an overflow. 

An example to this phenomenon is given in the 
appendix. 

Therefore, not to have the possibility of overflow, 
the numbers must be kept between -1/ > J~2 and 1/ • 

However, because of the implementation difficulty, 
correction is made to keep the numbers between -0,5 and 
0.5 (i.e., 600i 8 and 1777^.) 

CORRECTION TO DECUS NO. 8-143 

ORIGINAL CORRECTED 

x 1164 * 1164 

t?at. RTL 

CORRECTED 

5.2. Data Scaling 

All inputs X ± must be scaled such that 

X 1 4 1777 e’ Vi* 


n 




% +* 


- 2 - 

(if the length of the input data sequence is going to 
he at least doubled by extending it with zeros, input 
data must be scaled such that 4001^. 4*3777, 

appendix n 

An example to the error of the subroutine 
DECUS NO. 8-143 is given below. 


INPUT DATA (OCTAL) 

FOURIER ENERGY SPECTRUM OF THE INPUT DATA (DECIMAL) 

0000 

1411 

o************** ******** 

0100 

0256 

1**** <-ERROR (MUST BE 1017) 

0371 

0354- 

1 ****** 

1015 

0041 

1* 

1523 

00 

1 

2231 

0001 

# 

2656 

0000 

1 

3146 

0006 

1 <-ERROR (MUST BE 0000) 

3246 

0000 

1 

3146 

050 2 

1 ******** <- ERRQR (MUST BE jtytytyf) 

2656 

0000 

0 

2231 

0000 

1 

1523 

0000 

1 

1)015 

0000 

1 

J0371 

0000 

1 

0100 

0205 

.*** <— ERROR (MUST BE 0000) 

0000 

0000 

0000 

1 

0000 

SCALING FACTOR =2 T0004 

0000 

0000 

0000 

0000 

0000 

0000 

0000 

0000 

0000 

0000 

0000 

0000 



0000 


% 
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